//
// CUDA Sagan Simulator - a CUDA version fo the sagan Simulator Software
// version 1.0 - Alpha 01
// by TTYMan
//
// para compilar:  nvcc engine.cu -o teste1 -I/home/ttyman/NVIDIA_GPU_Computing_SDK/C/common/inc

#include <stdio.h>
#include <string.h>
#include <cutil.h>
//cutil.h ==> /home/ttyman/NVIDIA_GPU_Computing_SDK/C/common/inc 

#define BSIZE 256
#define BUFFER_MAX 2048      // buffer size for reading each line of initial conditions file.
#define offset      8        // each offset space is used to save one particle information:
                             //   0 - Pos X ==> array[index*offset+0]
                             //   1 - Pos Y ==> array[index*offset+1]
                             //   2 - Pos Z
                             //   3 - Mass
                             //   4 - Vel X
                             //   5 - Vel Y
                             //   6 - Vel Z
                             //   7 - not used
                             // this approuch is useful to take advantage of CUDA float3 and float4 datatypes.


float* LoadInitialConditionsData(char* FileName, int* N)
{
    FILE* fptr;

    if ((fptr = fopen(FileName, "r")))
    {
        char    LineBuffer[BUFFER_MAX];
        int     LineCount = 0;

        // The first line of the initial conditions file should contain the amount of particles information.
        // So, read this value and use it for allocate memory for the particles array.
        fgets (LineBuffer, BUFFER_MAX, fptr);
        sscanf(LineBuffer, "%d ", &LineCount);
        *N = LineCount; // 'N' value finally discovered!

        // allocate memory for the Particles array.
        float* ParticleArray = (float*)malloc( sizeof(float) * offset * LineCount );

        // read the file line by line, saving position and velocity data into arrays.
        // each file line represents one particle.
        for (int i=0; i<LineCount; i++)
        {
            fgets (LineBuffer, BUFFER_MAX, fptr);
            sscanf(LineBuffer, "%f %f %f %f %f %f %f " , &ParticleArray[i*offset+0],      /* PosX     */
                                                         &ParticleArray[i*offset+1],      /* PosY     */
                                                         &ParticleArray[i*offset+2],      /* PosZ     */
                                                         &ParticleArray[i*offset+3],      /* Mass     */
                                                         &ParticleArray[i*offset+4],      /* VelX     */
                                                         &ParticleArray[i*offset+5],      /* VelY     */
                                                         &ParticleArray[i*offset+6]);     /* VelZ     */
                                                          ParticleArray[i*offset+7]=0.0f; /* Not Used */

        }

        //close the file and release ay used resource
        fclose(fptr);

        return ParticleArray;

    }
    else
    {
        //Error! Couldn't open the file!
        printf("\nError opening file %s!\n", FileName);
        return 0;
    }
}

__device__ float3
ParticleParticleInteraction(float4 ParticleI, float4 ParticleJ, float3 AccelerationI, float SofteningFactor)
{
	
    //float Force;
    float XDist;
    float YDist;
    float ZDist;
    float XT;
    float YT;
    float ZT;
    float DistanceSQR;
    float Rinv;
    float Rinv3;
 
    //float eps = 4.0f/N; // softening factor
	// SofteningFactor = eps
	
	XDist = ParticleI.x - ParticleJ.x;
    XT = XDist * XDist;

    YDist = ParticleI.y - ParticleJ.y;
    YT = YDist * YDist;

    ZDist = ParticleI.z - ParticleJ.z;
    ZT = ZDist * ZDist;

    DistanceSQR = XT + YT + ZT;

    Rinv = 1.0f / sqrtf( DistanceSQR + (SofteningFactor * SofteningFactor) );

    Rinv3 = cbrtf(Rinv);

    AccelerationI.x -= ParticleJ.w * XDist * Rinv3;
	AccelerationI.y -= ParticleJ.w * YDist * Rinv3;
	AccelerationI.z -= ParticleJ.w * ZDist * Rinv3;

    return AccelerationI; 
}

__device__ float3
tile_calculation(float4 myPosition, float3 Acceleration)
{
	extern __shared__ float4 shPosition[];
	
	#pragma unroll 8
	for (unsigned int i = 0; i < BSIZE; i++)
		Acceleration = ParticleParticleInteraction(myPosition, shPosition[i], Acceleration, 0.0f);
		
	return Acceleration;
}

__global__ void 
galaxyKernel(float4* pos, float4 * pdata, unsigned int width, 
			 unsigned int height, float step, int apprx, int tile_offset)
{
	// shared memory
	extern __shared__ float4 shPosition[];
	
	// index of my body	
	unsigned int x = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int y = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned int pLoc = y * width + x;
    unsigned int vLoc = width * height + pLoc;
	
    // starting index of position array
    unsigned int start = ( (width * height) / apprx ) * tile_offset;
	
	float4 myPosition = pdata[pLoc];
	float4 myVelocity = pdata[vLoc];

	float3 Acceleration = {0.0f, 0.0f, 0.0f};

	unsigned int idx = 0;
	unsigned int loop = ((width * height) / apprx ) / BSIZE;
	for (int i = 0; i < loop; i++)
	{
		idx = threadIdx.y * blockDim.x + threadIdx.x;
		shPosition[idx] = pdata[idx + start + BSIZE * i];

		__syncthreads();
		
		Acceleration = tile_calculation(myPosition, Acceleration);
		
		__syncthreads();		
	}
    	
    // update velocity with above acc
    myVelocity.x += Acceleration.x * step;// * 2.0f;
    myVelocity.y += Acceleration.y * step;// * 2.0f;
    myVelocity.z += Acceleration.z * step;// * 2.0f;
    
    //myVelocity.x *= damping;
    //myVelocity.y *= damping;
    //myVelocity.z *= damping;
    
    // update position
    myPosition.x += myVelocity.x * step;
    myPosition.y += myVelocity.y * step;
    myPosition.z += myVelocity.z * step;
        
    __syncthreads();
    
    // update device memory
	pdata[pLoc] = myPosition;
	pdata[vLoc] = myVelocity;
    
	// update vbo
	pos[pLoc] = make_float4(myPosition.x, myPosition.y, myPosition.z, 1.0f);
	pos[vLoc] = myVelocity;
}

// void CallCUDAKernel(float4* pos, float4 * pdata, int width, int height, float step, int apprx, int tile_offset, int N)
// {
    // dim3 block(16, 16, 1);
    // dim3 grid(width / block.x, height / block.y, 1);
    // int sharedMemSize = BSIZE * sizeof(float4);
	
	// float4* d_pos;
	// float4* d_pdata;
	// int    d_width;
	// int    d_height;
	// float  d_step;
	// float  d_apprx;
	// int    d_tile_offset;
	
	//allocate device memory
	// cudaMalloc( (void**)&d_pos,          sizeof(float4)   );
	// cudaMalloc( (void**)&d_pdata,        sizeof(float4)*N );
	// cudaMalloc( (void**)&d_width,        sizeof(int)      );
	// cudaMalloc( (void**)&d_height,       sizeof(int)      );
	// cudaMalloc( (void**)&d_step,         sizeof(float)    );
	// cudaMalloc( (void**)&d_apprx,        sizeof(float)    );
	// cudaMalloc( (void**)&d_tile_offset,  sizeof(int)      );
	
	//send its copies to the device memory
	// cudaMemcpy( d_pos,         pos,         sizeof(float4),   cudaMemcpyHostToDevice );
	// cudaMemcpy( d_pdata,       pdata,       sizeof(float4)*N, cudaMemcpyHostToDevice );
	// cudaMemcpy( d_width,       width,       sizeof(int),      cudaMemcpyHostToDevice );
	// cudaMemcpy( d_height,      height,      sizeof(int),      cudaMemcpyHostToDevice );
	// cudaMemcpy( d_step,        step,        sizeof(float),    cudaMemcpyHostToDevice );
	// cudaMemcpy( d_apprx,       apprx,       sizeof(float),    cudaMemcpyHostToDevice );
	// cudaMemcpy( d_tile_offset, tile_offset, sizeof(int),      cudaMemcpyHostToDevice );
	
	
    //galaxyKernel<<<grid, block, sharedMemSize>>>(pos, pdata, width, height, step, apprx, tile_offset);
// }



int main(int argc, char **argv)
{
    printf("\nOla, mundo!!!\n");

    char    Arquivo[] = "./data/esfera395CUDA.ini";
    int     N = 0;
    float*  ParticleArray;
	float4* d_ParticleArray;
	float4* headPtr;
	float4* d_headPtr;
	unsigned int       width = 32;
	unsigned int    d_width;
	unsigned int       height = 32;
	unsigned int    d_height;
	float     step;
	float*  d_step;
	//float apprx, d_apprx;
	int      tile_offset;
	int*   d_tile_offset;

    dim3 block(16, 16, 1);
    dim3 grid(width / block.x, height / block.y, 1);
    int sharedMemSize = BSIZE * sizeof(float4);
	
	
    ParticleArray = LoadInitialConditionsData(Arquivo, &N);

	cudaMalloc( (void**)&d_headPtr,         sizeof(float4)   );
	cudaMalloc( (void**)&d_ParticleArray, sizeof(float4)*N );
	cudaMalloc( (void**)&d_width,         sizeof(int)      );
	cudaMalloc( (void**)&d_height,        sizeof(int)      );
	cudaMalloc( (void**)&d_step,          sizeof(float)    );
	//cudaMalloc( (void**)&d_apprx,         sizeof(float)    );
	cudaMalloc( (void**)&d_tile_offset,   sizeof(int)      );
	
	cudaMemcpy( d_ParticleArray, &ParticleArray, sizeof(float4)*N, cudaMemcpyHostToDevice );
	cudaMemcpy( d_width,         &width,         sizeof(unsigned int),      cudaMemcpyHostToDevice );
	cudaMemcpy( d_height,        &height,        sizeof(unsigned int),      cudaMemcpyHostToDevice );
	cudaMemcpy( d_step,          &step,          sizeof(float),    cudaMemcpyHostToDevice );
	cudaMemcpy( d_tile_offset,   &tile_offset,   sizeof(int),      cudaMemcpyHostToDevice );
	
	
    //                                            width, height, time step, approximation, tile_offset, N"amount of particles"
	//CallCUDAKernel( dptr, (float4*)ParticleArray, 256, N / 256, 0.0001f, 1, 16, N);	
	
	galaxyKernel<<<grid, block, sharedMemSize>>>(d_headPtr, 
	d_ParticleArray, 
	d_width, 
	d_height, 
	d_step, 1, d_tile_offset);







	// teste teste
    int devID;
    cudaDeviceProp props;
    CUDA_SAFE_CALL( cudaGetDevice(&devID) );
    CUDA_SAFE_CALL( cudaGetDeviceProperties(&props, devID) );


    printf("\nValor de N: %d\n", N);

    // // teste teste
    // for(int i=0; i<N; i++)
            // printf("main i=%d: %f %f %f %f %f %f %f \n" , i,          ParticleArray[i*offset+0],      /* PosX     */
                                                         // ParticleArray[i*offset+1],      /* PosY     */
                                                         // ParticleArray[i*offset+2],      /* PosZ     */
                                                         // ParticleArray[i*offset+3],      /* Mass     */
                                                         // ParticleArray[i*offset+4],      /* VelX     */
                                                         // ParticleArray[i*offset+5],      /* VelY     */
                                                         // ParticleArray[i*offset+6]);     /* VelZ     */

    printf("\n\nYou ran me by typing: %s\n", argv[0]);
	printf("\nPlaca: %s\n", props.name);
    return (0);
}
